Consider the cohort as a single group to start. Do not stratify by case-control status.

Stack the heatmaps for CD4. Are there obvious patterns to the antigens?
Stack the heatmaps for NOTCD4. Are there patterns? Are these patterns different from CD4?

Calculate FS/PFS for each stimulation in CD4 and compare side-by-side in boxplots and test by ANOVA. Is there a clear winner? Calculate FS/PFS for each stimulation in NOTCD4 and compare with boxplot and ANOVA. Is there a winner? Is this different than CD4?

It may also make sense to make a corrplot of the FS/PFS from the different stims. This is asking a different question though, i.e. whether a patient’s cytokine responses tend to be consistently elevated/not-elevated across multiple stims.

Regress FS/PFS against age (linear). Is there an effect of age on T cell functionality? Regress FS/PFS against days from symptom onset (linear). Is there a decrease in functionality with time? Stratify FS/PFS by sex and test by t-test (or Mann-Whitney, though sample size is large enough). Is there a difference?

Only then would you start asking the case-control questions.

library(here)
library(tidyverse)
library(COMPASS)
library(ggplot2)
library(data.table)
library(broom)
library(psych)
library(corrplot)
library(ggpubr)
library(knitr)
library(kableExtra)
library(cowplot)
source(here::here("scripts/my_pheatmap_and_plotMeanGamma.R"))
source(here::here("scripts/20200604_Helper_Functions.R"))
save_output <- FALSE

Read in data

# Next 4 lines from Run_COMPASS.R
stims_for_compass_runs <- c("VEMP", "Spike 1", "Spike 2", "NCAP")
parent_nodes_for_compass_runs <- c("4+", "NOT4+")
stims_for_compass_runs_rep <- rep(stims_for_compass_runs, each = length(parent_nodes_for_compass_runs))
parent_nodes_for_compass_runs_rep <- rep(parent_nodes_for_compass_runs, times = length(stims_for_compass_runs))

crList <- purrr::pmap(.l = list(parent_nodes_for_compass_runs_rep,
                                stims_for_compass_runs_rep),
                      .f = function(parent, currentStim) {
                        currentStimForFile <- gsub(" ", "_", currentStim)
                        crPath <- here::here(sprintf("out/CompassOutput/%s/%s/COMPASSResult_%s_%s.rds", parent, currentStimForFile, parent, currentStimForFile))
                        readRDS(crPath)
                      }) %>% 
  set_names(gsub("+", "", gsub(" ", "_", paste0(parent_nodes_for_compass_runs_rep, "_", stims_for_compass_runs_rep)), fixed=TRUE))
names(crList)
## [1] "4_VEMP"       "NOT4_VEMP"    "4_Spike_1"    "NOT4_Spike_1" "4_Spike_2"   
## [6] "NOT4_Spike_2" "4_NCAP"       "NOT4_NCAP"
CD4RunNames <- c("4_Spike_1", "4_Spike_2", "4_NCAP", "4_VEMP")
NotCD4RunNames <- c("NOT4_Spike_1", "NOT4_Spike_2", "NOT4_NCAP", "NOT4_VEMP")

Stacked heatmaps

plotCOMPASSResultStack(crList[CD4RunNames],
                       row_annotation = c("Stim"),
                       variable = "Stim",
                       show_rownames = FALSE,
                       main = "CD4 COMPASS Runs",
                       fontsize = 15, fontsize_row = 15, fontsize_col = 13)
## The 'threshold' filter has removed 6 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&IL17a Ax700&!IL2 PE&!IL4/5/13 APC&!TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&!TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&IL4/5/13 APC&!TNFa FITC

plotCOMPASSResultStack(crList[NotCD4RunNames],
                       row_annotation = c("Stim"),
                       variable = "Stim",
                       show_rownames = FALSE,
                       main = "Not-CD4 COMPASS Runs",
                       fontsize = 15, fontsize_row = 15, fontsize_col = 13)
## The 'threshold' filter has removed 2 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&IL17a Ax700&!IL2 PE&IL4/5/13 APC&!TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC

plotCOMPASSResultStack(crList[c(CD4RunNames, NotCD4RunNames)],
                       row_annotation = c("Stim"),
                       variable = "Stim",
                       show_rownames = FALSE,
                       main = "CD4 and Not-CD4 COMPASS Runs",
                       fontsize = 15, fontsize_row = 15, fontsize_col = 13)
## The 'threshold' filter has removed 3 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&IL4/5/13 APC&!TNFa FITC

if(save_output) {
  png(filename=here::here("out/PostCompassPlots/20200613_CD4_COMPASS_Heatmap_All_Stims.png"),
      width=700, height=1000)
  plotCOMPASSResultStack(crList[CD4RunNames],
                         row_annotation = c("Stim"),
                         variable = "Stim",
                         show_rownames = FALSE,
                         main = "CD4 COMPASS Runs",
                         fontsize = 14, fontsize_row = 13, fontsize_col = 13)
  dev.off()
  
  png(filename=here::here("out/PostCompassPlots/20200613_NotCD4_COMPASS_Heatmap_All_Stims.png"),
      width=700, height=1000)
  plotCOMPASSResultStack(crList[NotCD4RunNames],
                         row_annotation = c("Stim"),
                         variable = "Stim",
                         show_rownames = FALSE,
                         main = "Not-CD4 COMPASS Runs",
                         fontsize = 14, fontsize_row = 13, fontsize_col = 13)
  dev.off()
  
  png(filename=here::here("out/PostCompassPlots/20200613_CD4_and_NotCD4_COMPASS_Heatmap_All_Stims.png"),
      width=900, height=2000)
  plotCOMPASSResultStack(crList[c(CD4RunNames, NotCD4RunNames)],
                         row_annotation = c("Stim"),
                         variable = "Stim",
                         show_rownames = FALSE,
                         main = "CD4 and Not-CD4 COMPASS Runs",
                         fontsize = 14, fontsize_row = 13, fontsize_col = 13)
  dev.off()
}

FS/PFS across Stims

fs_pfs_df <- lapply(c(CD4RunNames, NotCD4RunNames), function(n) {
  as.data.frame(COMPASS::FunctionalityScore(crList[[n]])) %>%
    rownames_to_column("SAMPLE ID") %>% 
    rename(FS = `COMPASS::FunctionalityScore(crList[[n]])`) %>% 
    left_join(as.data.frame(COMPASS::PolyfunctionalityScore(crList[[n]])) %>%
      rownames_to_column("SAMPLE ID") %>% 
      rename(PFS = `COMPASS::PolyfunctionalityScore(crList[[n]])`),
      by = "SAMPLE ID") %>% 
    mutate(COMPASS_run = n)
  } ) %>% 
  data.table::rbindlist() %>% 
  separate(COMPASS_run, into = c("parent", "STIM"), remove = F, extra = "merge") %>% 
  mutate(STIM = factor(STIM, levels = c("Spike_1", "Spike_2", "NCAP", "VEMP")))#,
         #`SAMPLE ID` = factor(`SAMPLE ID`, levels = unique(`SAMPLE ID`)))

table(fs_pfs_df$STIM, fs_pfs_df$`SAMPLE ID`) # remove 07B for quade.test in order to have complete block design
##          
##           07B 102C 103C 106C 109C 112C 113C 117C 121C 122C 127C 128C 12C 13
##   Spike_1   2    2    2    2    2    2    2    2    2    2    2    2   2  2
##   Spike_2   0    2    2    2    2    2    2    2    2    2    2    2   2  2
##   NCAP      2    2    2    2    2    2    2    2    2    2    2    2   2  2
##   VEMP      2    2    2    2    2    2    2    2    2    2    2    2   2  2
##          
##           130C 131C 133C 136C 139C 141C 142C 143C 147C 148C 150C 15514 15518
##   Spike_1    2    2    2    2    2    2    2    2    2    2    2     2     2
##   Spike_2    2    2    2    2    2    2    2    2    2    2    2     2     2
##   NCAP       2    2    2    2    2    2    2    2    2    2    2     2     2
##   VEMP       2    2    2    2    2    2    2    2    2    2    2     2     2
##          
##           15529 15530 15545 15548 15554 15575 157C 159C 161C 162C 164C 169C
##   Spike_1     2     2     2     2     2     2    2    2    2    2    2    2
##   Spike_2     2     2     2     2     2     2    2    2    2    2    2    2
##   NCAP        2     2     2     2     2     2    2    2    2    2    2    2
##   VEMP        2     2     2     2     2     2    2    2    2    2    2    2
##          
##           170C 1C 21C 23 25 25C 29C 36C 38C 51C 56C 59C 69C 6C 76C 80C 90C 93C
##   Spike_1    2  2   2  2  2   2   2   2   2   2   2   2   2  2   2   2   2   2
##   Spike_2    2  2   2  2  2   2   2   2   2   2   2   2   2  2   2   2   2   2
##   NCAP       2  2   2  2  2   2   2   2   2   2   2   2   2  2   2   2   2   2
##   VEMP       2  2   2  2  2   2   2   2   2   2   2   2   2  2   2   2   2   2
##          
##           96C BWT20 CLT1 HS09 HS10 HS8
##   Spike_1   2     2    2    2    2   2
##   Spike_2   2     2    2    2    2   2
##   NCAP      2     2    2    2    2   2
##   VEMP      2     2    2    2    2   2
fs_pfs_df.complete <- fs_pfs_df %>%
             dplyr::filter(`SAMPLE ID` != "07B") %>% 
             mutate(`SAMPLE ID` = factor(`SAMPLE ID`, levels = unique(`SAMPLE ID`))) %>% 
  arrange(STIM, `SAMPLE ID`)
fs_pfs_df.complete.4 <- fs_pfs_df.complete %>% dplyr::filter(parent == "4")
fs_pfs_df.complete.not4 <- fs_pfs_df.complete %>% dplyr::filter(parent == "NOT4")

# Perform statistical tests for Functionality Score across Stims
# Quade test is to wilcoxon signed rank test like ANOVA is to t-test
quade_fs_4 <- quade.test(FS ~ STIM | `SAMPLE ID`, data = fs_pfs_df.complete.4)
quade_fs_4
## 
##  Quade test
## 
## data:  FS and STIM and SAMPLE ID
## Quade F = 73.3, num df = 3, denom df = 183, p-value < 2.2e-16
quade_fs_not4 <- quade.test(FS ~ STIM | `SAMPLE ID`, data = fs_pfs_df.complete.not4)
quade_fs_not4
## 
##  Quade test
## 
## data:  FS and STIM and SAMPLE ID
## Quade F = 41.83, num df = 3, denom df = 183, p-value < 2.2e-16
pw_fs_4 <- pairwise.wilcox.test(fs_pfs_df.complete.4$FS,
                      fs_pfs_df.complete.4$STIM,
                      p.adjust.method="bonferroni",
                      paired=TRUE)
pw_fs_4$p.value
##              Spike_1      Spike_2         NCAP
## Spike_2 1.632616e-01           NA           NA
## NCAP    4.082570e-02 3.720152e-06           NA
## VEMP    6.244948e-11 4.659364e-11 6.244948e-11
broom::tidy(pw_fs_4)
## # A tibble: 6 x 3
##   group1  group2   p.value
##   <chr>   <chr>      <dbl>
## 1 Spike_2 Spike_1 1.63e- 1
## 2 NCAP    Spike_1 4.08e- 2
## 3 VEMP    Spike_1 6.24e-11
## 4 NCAP    Spike_2 3.72e- 6
## 5 VEMP    Spike_2 4.66e-11
## 6 VEMP    NCAP    6.24e-11
pw_fs_not4 <- pairwise.wilcox.test(fs_pfs_df.complete.not4$FS,
                      fs_pfs_df.complete.not4$STIM,
                      p.adjust.method="bonferroni",
                      paired=TRUE)
pw_fs_not4$p.value
##              Spike_1      Spike_2         NCAP
## Spike_2 1.000000e+00           NA           NA
## NCAP    7.599030e-09 1.461474e-09           NA
## VEMP    1.984327e-02 2.881666e-01 4.411433e-10
broom::tidy(pw_fs_not4)
## # A tibble: 6 x 3
##   group1  group2   p.value
##   <chr>   <chr>      <dbl>
## 1 Spike_2 Spike_1 1.00e+ 0
## 2 NCAP    Spike_1 7.60e- 9
## 3 VEMP    Spike_1 1.98e- 2
## 4 NCAP    Spike_2 1.46e- 9
## 5 VEMP    Spike_2 2.88e- 1
## 6 VEMP    NCAP    4.41e-10
# Perform statistical tests for PolyFunctionality Score across Stims
quade_pfs_4 <- quade.test(PFS ~ STIM | `SAMPLE ID`, data = fs_pfs_df.complete.4)
quade_pfs_4
## 
##  Quade test
## 
## data:  PFS and STIM and SAMPLE ID
## Quade F = 63.947, num df = 3, denom df = 183, p-value < 2.2e-16
quade_pfs_not4 <- quade.test(PFS ~ STIM | `SAMPLE ID`, data = fs_pfs_df.complete.not4)
quade_pfs_not4
## 
##  Quade test
## 
## data:  PFS and STIM and SAMPLE ID
## Quade F = 32.722, num df = 3, denom df = 183, p-value < 2.2e-16
pw_pfs_4 <- pairwise.wilcox.test(fs_pfs_df.complete.4$PFS,
                      fs_pfs_df.complete.4$STIM,
                      p.adjust.method="bonferroni",
                      paired=TRUE)
pw_pfs_4$p.value
##              Spike_1      Spike_2         NCAP
## Spike_2 3.831211e-02           NA           NA
## NCAP    1.000000e+00 6.717394e-04           NA
## VEMP    7.960767e-11 4.659364e-11 6.556256e-11
# broom::tidy(pw_pfs_4)

pw_pfs_not4 <- pairwise.wilcox.test(fs_pfs_df.complete.not4$PFS,
                      fs_pfs_df.complete.not4$STIM,
                      p.adjust.method="bonferroni",
                      paired=TRUE)
pw_pfs_not4$p.value
##              Spike_1      Spike_2         NCAP
## Spike_2 1.000000e+00           NA           NA
## NCAP    2.137504e-08 1.518882e-08           NA
## VEMP    1.000000e+00 1.000000e+00 3.142649e-09
# broom::tidy(pw_pfs_not4)

# Now plot the results
# Note that the dataset used in statistics excludes 07B, but 07B is plotted below (and the median bar is calculated including 07B)

parent.labs <- c("CD4", "Not-CD4")
names(parent.labs) <- c("4", "NOT4")

plot_fs_pfs_vs_stim <- function(fs_or_pfs) {
  ggplot(fs_pfs_df,
         aes(STIM, !!as.name(fs_or_pfs))) +
    theme_bw(base_size = 22) +
    geom_line(aes(group = `SAMPLE ID`), color="grey") +
    geom_point(color="grey") +
    facet_grid(. ~ parent,
               labeller = labeller(parent = parent.labs)) +
    geom_errorbarh(data = fs_pfs_df %>% group_by(parent, STIM) %>% summarise(med = median(!!as.name(fs_or_pfs))),
                   aes(y = med,
                       xmax = 1.5 + 0.35,
                       xmin = 1.5 - 0.35,
                       height = 0),
                   position=position_dodge(width=0), color = "black") +
    labs(title = sprintf("%s vs Stim", fs_or_pfs),
         y = fs_or_pfs) +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_text(size=22),
          axis.text.y = element_text(color="black", size=22),
          axis.text.x = element_text(color="black", size=22),
          plot.title = element_text(hjust = 0.5),
          panel.grid = element_blank(),
          legend.position = "none",
          plot.margin = margin(1.3, 0.2, 0, 0.2, "cm")) +
    scale_x_discrete(labels=c(c("Spike_1" = "Spike 1", "Spike_2" = "Spike 2")),
                     expand = c(0.1,0.1))
}

plot_fs_pfs_vs_stim("FS")

# Remembering these p-values are adjusted
as.data.frame(pw_fs_4$p.value) %>% 
  rownames_to_column('tmp') %>% 
  mutate_all(~cell_spec(.x, color = ifelse(!is.na(.x) & .x < 0.05, "red", "black"))) %>%
  column_to_rownames('tmp') %>% 
  kable(format = "html", escape = F, row.names = T) %>%
  kable_styling("striped", full_width = F)
Spike_1 Spike_2 NCAP
Spike_2 0.163261630177697 NA NA
NCAP 0.0408257009089326 3.72015208085264e-06 NA
VEMP 6.24494803673655e-11 4.65936369334732e-11 6.24494803673655e-11
as.data.frame(pw_fs_not4$p.value) %>% 
  rownames_to_column('tmp') %>% 
  mutate_all(~cell_spec(.x, color = ifelse(!is.na(.x) & .x < 0.05, "red", "black"))) %>%
  column_to_rownames('tmp') %>% 
  kable(format = "html", escape = F, row.names = T) %>%
  kable_styling("striped", full_width = F)
Spike_1 Spike_2 NCAP
Spike_2 1 NA NA
NCAP 7.5990300522029e-09 1.46147350795162e-09 NA
VEMP 0.0198432725984609 0.288166634062905 4.41143306237813e-10
plot_fs_pfs_vs_stim("PFS")

as.data.frame(pw_pfs_4$p.value) %>% 
  rownames_to_column('tmp') %>% 
  mutate_all(~cell_spec(.x, color = ifelse(!is.na(.x) & .x < 0.05, "red", "black"))) %>%
  column_to_rownames('tmp') %>% 
  kable(format = "html", escape = F, row.names = T) %>%
  kable_styling("striped", full_width = F)
Spike_1 Spike_2 NCAP
Spike_2 0.0383121094257177 NA NA
NCAP 1 0.000671739390258177 NA
VEMP 7.96076669051054e-11 4.65936369334732e-11 6.55625616466633e-11
as.data.frame(pw_pfs_not4$p.value) %>% 
  rownames_to_column('tmp') %>% 
  mutate_all(~cell_spec(.x, color = ifelse(!is.na(.x) & .x < 0.05, "red", "black"))) %>%
  column_to_rownames('tmp') %>% 
  kable(format = "html", escape = F, row.names = T) %>%
  kable_styling("striped", full_width = F)
Spike_1 Spike_2 NCAP
Spike_2 1 NA NA
NCAP 2.13750379422692e-08 1.51888161479098e-08 NA
VEMP 1 1 3.14264895925017e-09

As we saw in the heatmaps, NAP and Spike 1 consistently have the highest FS/PFS scores, and VEMP has the lowest.

How correlated are FS/PFS across stims?

fs_pfs_df_wide <- fs_pfs_df %>%  
               pivot_wider(id_cols = c(`SAMPLE ID`),
                           names_from = c(STIM, parent),
                           values_from = c(FS, PFS))

pairs.panels(fs_pfs_df_wide %>% 
               dplyr::select(-c(`SAMPLE ID`)) %>% 
               dplyr::select(matches("^FS.*_4")), 
             method = "pearson", # correlation method, "r" displayed in upper off diagonals. Not r^2.
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots,
             lm = T,
             main = "CD4 FS Correlation across Stims")

pairs.panels(fs_pfs_df_wide %>% 
               dplyr::select(-c(`SAMPLE ID`)) %>% 
               dplyr::select(matches("^FS.*_NOT4")), 
             method = "pearson", # correlation method, "r" displayed in upper off diagonals. Not r^2.
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots,
             lm = T,
             main = "Not-CD4 FS Correlation across Stims")

pairs.panels(fs_pfs_df_wide %>% 
               dplyr::select(-c(`SAMPLE ID`)) %>% 
               dplyr::select(matches("^PFS.*_4")), 
             method = "pearson", # correlation method, "r" displayed in upper off diagonals. Not r^2.
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots,
             lm = T,
             main = "CD4 FS Correlation across Stims")

pairs.panels(fs_pfs_df_wide %>% 
               dplyr::select(-c(`SAMPLE ID`)) %>% 
               dplyr::select(matches("^PFS.*_NOT4")), 
             method = "pearson", # correlation method, "r" displayed in upper off diagonals. Not r^2.
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots,
             lm = T,
             main = "Not-CD4 PFS Correlation across Stims")

M <- cor(fs_pfs_df_wide %>% 
               dplyr::select(-c(`SAMPLE ID`)) %>% 
               dplyr::select(matches("PFS")),
         use = "complete.obs") # missing values are handled by casewise deletion
# mat : is a matrix of data
# ... : further arguments to pass to the native R cor.test function
cor.mtest <- function(mat, ...) {
    mat <- as.matrix(mat)
    n <- ncol(mat)
    p.mat<- matrix(NA, n, n)
    diag(p.mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            tmp <- cor.test(mat[, i], mat[, j], ...)
            p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
        }
    }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}
getOption("na.action") # Default behavior of cor.test
## [1] "na.omit"
# matrix of the p-value of the correlation
p.mat <- cor.mtest(fs_pfs_df_wide %>% 
               dplyr::select(-c(`SAMPLE ID`)) %>% 
               dplyr::select(matches("PFS")))
head(p.mat[, 1:5])
##                  PFS_Spike_1_4 PFS_Spike_2_4   PFS_NCAP_4   PFS_VEMP_4
## PFS_Spike_1_4     0.000000e+00  3.557519e-07 1.252908e-08 0.0003300951
## PFS_Spike_2_4     3.557519e-07  0.000000e+00 1.010421e-05 0.0196248181
## PFS_NCAP_4        1.252908e-08  1.010421e-05 0.000000e+00 0.0018545939
## PFS_VEMP_4        3.300951e-04  1.962482e-02 1.854594e-03 0.0000000000
## PFS_Spike_1_NOT4  2.113868e-09  4.275608e-03 6.103202e-03 0.0173764534
## PFS_Spike_2_NOT4  2.010616e-01  5.250470e-04 2.707909e-01 0.1543216526
##                  PFS_Spike_1_NOT4
## PFS_Spike_1_4        2.113868e-09
## PFS_Spike_2_4        4.275608e-03
## PFS_NCAP_4           6.103202e-03
## PFS_VEMP_4           1.737645e-02
## PFS_Spike_1_NOT4     0.000000e+00
## PFS_Spike_2_NOT4     2.657072e-01
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(M,
         method="color", col=col(200),  
         type="upper", order="hclust", 
         addCoef.col = "black", # Add coefficient of correlation
         tl.col="black", tl.srt=45, #Text label color and rotation
         # Combine with significance
         p.mat = p.mat, sig.level = 0.01, insig = "blank", 
         # hide correlation coefficient on the principal diagonal
         diag=FALSE)

# significant coefficients are colored

Observations:
(1) Note that FS and PFS scores for the same condition tend to be highly correlated (not shown). Not surprising.
(2) Within PFS only, there is a mixed bag of significant associations occuring “between CD4 and Not-CD4 for the same condition” and “between different stims”. Tends to be much lower r than those for (1) though.

In general, we can’t rely on a particular condition to predict the FS/PFS of all the other conditions. Each STIM and co-receptor subset provides different information. (If the former were true, it might have simplified our upcoming analyses)

Regress FS/PFS vs Demographics

# Read in the final patient manifest with more complete data  
data_manifest <- read.csv(here::here("data/Seshadri_HAARVI_PBMC_manifest_merged_11June2020.csv"), check.names = F, stringsAsFactors = F)

setdiff(sort(unique(fs_pfs_df$`SAMPLE ID`)), sort(unique(data_manifest$`Record ID`)))
## [1] "07B"   "BWT20" "CLT1"  "HS09"
setdiff(sort(unique(data_manifest$`Record ID`)), sort(unique(fs_pfs_df$`SAMPLE ID`)))
## [1] "116C"   "37C"    "551432" "75C"    "HS9"
# Remember 116C and 37C didn't get run through COMPASS
# 75C was not run at all
# 07B/551432 didn't get run for Spike 2
fs_pfs_df_with_manifest <- fs_pfs_df %>%
  mutate("Record ID" = dplyr::recode(`SAMPLE ID`,
                                     "07B" = "551432",
                                     "HS09" = "HS9")) %>% # also, 75C did not get run
  full_join(data_manifest, by = c("Record ID" = "Record ID")) %>% 
  dplyr::filter(!(`Record ID` %in% c("116C", "37C", "75C"))) %>% 
  mutate(Days_Symptom_Onset_to_Visit_1 = as.numeric(`Days symptom onset to visit 1`),
         Cohort = factor(ifelse(Cohort %in%
                                                c(NA, "Healthy control", "Healthy control 2017-2018"), "Healthy", Cohort),
                                              levels = c("Healthy", "Non-hospitalized", "Hospitalized")))
## Warning: NAs introduced by coercion
fs_pfs_df_with_manifest %>% dplyr::filter(is.na(FS) | is.na(`Record ID`) | is.na(`SAMPLE ID`) | is.na(STIM) | is.na(Cohort))
##  [1] SAMPLE ID                     FS                           
##  [3] PFS                           COMPASS_run                  
##  [5] parent                        STIM                         
##  [7] Record ID                     Sample ID                    
##  [9] Collection date               Cell count                   
## [11] Cohort                        Age                          
## [13] Sex                           Race                         
## [15] Hispanic?                     Days symptom onset to visit 1
## [17] Pair ID                       Race_v2                      
## [19] Days_Symptom_Onset_to_Visit_1
## <0 rows> (or 0-length row.names)

FS/PFS vs Age

# Note that p-values on ggscatter plots below are *not* adjusted, e.g.:
broom::tidy(lm(FS ~ Age, data = fs_pfs_df_with_manifest %>% dplyr::filter(parent == "4" & STIM == "Spike_2")))
## # A tibble: 2 x 5
##   term        estimate std.error statistic     p.value
##   <chr>          <dbl>     <dbl>     <dbl>       <dbl>
## 1 (Intercept) 0.0397    0.00704       5.64 0.000000529
## 2 Age         0.000230  0.000128      1.80 0.0772
ggscatter(fs_pfs_df_with_manifest,
          x = "Age", y = "FS",
          color = "STIM", palette = "jco",
          add = "reg.line") +
  facet_grid(parent ~ STIM) +
  stat_cor()
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 16 rows containing non-finite values (stat_smooth).
## Warning: Removed 16 rows containing non-finite values (stat_cor).
## Warning: Removed 16 rows containing missing values (geom_point).

ggscatter(fs_pfs_df_with_manifest,
          x = "Age", y = "PFS",
          color = "STIM", palette = "jco",
          add = "reg.line") +
  facet_grid(parent ~ STIM) +
  stat_cor()
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 16 rows containing non-finite values (stat_smooth).
## Warning: Removed 16 rows containing non-finite values (stat_cor).
## Warning: Removed 16 rows containing missing values (geom_point).

FS/PFS vs Days from Symptom Onset to Visit

# Note that p-values on ggscatter plots below are *not* adjusted, e.g.:
broom::tidy(lm(FS ~ Days_Symptom_Onset_to_Visit_1, data = fs_pfs_df_with_manifest %>% dplyr::filter(parent == "4" & STIM == "Spike_2")))
## # A tibble: 2 x 5
##   term                          estimate std.error statistic  p.value
##   <chr>                            <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)                   0.0469    0.00600       7.83 1.68e-10
## 2 Days_Symptom_Onset_to_Visit_1 0.000129  0.000114      1.13 2.63e- 1
# Note that the Healthies don't have values for Days_Symptom_Onset_to_Visit_1 and will get filtered out below
ggscatter(fs_pfs_df_with_manifest,
          x = "Days_Symptom_Onset_to_Visit_1", y = "FS",
          color = "STIM", palette = "jco",
          add = "reg.line") +
  facet_grid(parent ~ STIM) +
  stat_cor()
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 46 rows containing non-finite values (stat_smooth).
## Warning: Removed 46 rows containing non-finite values (stat_cor).
## Warning: Removed 46 rows containing missing values (geom_point).

ggscatter(fs_pfs_df_with_manifest,
          x = "Days_Symptom_Onset_to_Visit_1", y = "PFS",
          color = "STIM", palette = "jco",
          add = "reg.line") +
  facet_grid(parent ~ STIM) +
  stat_cor()
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 46 rows containing non-finite values (stat_smooth).
## Warning: Removed 46 rows containing non-finite values (stat_cor).
## Warning: Removed 46 rows containing missing values (geom_point).

For the CD4 FS vs Age plots, there is a positive association of Age with FS for Spike 1 (at least).

FS/PFS vs Sex

# Note that p-values on ggboxplot plots below are *not* adjusted, e.g.:
broom::tidy(wilcox.test(FS ~ Sex, data = fs_pfs_df_with_manifest %>% dplyr::filter(parent == "4" & STIM == "Spike_2")))
## # A tibble: 1 x 4
##   statistic p.value method                 alternative
##       <dbl>   <dbl> <chr>                  <chr>      
## 1       369   0.239 Wilcoxon rank sum test two.sided
# Note that future ggpubr update should make outlier.shape=NA occur by default when jitter is added, but include for now
ggboxplot(fs_pfs_df_with_manifest %>% 
            dplyr::filter(!is.na(Sex)),
          x = "Sex", y = "FS",
          color = "Sex", palette = "jco",
          add = "jitter",
          outlier.shape = NA) + 
  facet_grid(parent ~ STIM) +
  stat_compare_means(aes(group = Sex))

ggboxplot(fs_pfs_df_with_manifest %>% 
            dplyr::filter(!is.na(Sex)),
          x = "Sex", y = "PFS",
          color = "Sex", palette = "jco",
          add = "jitter",
          outlier.shape = NA) + 
  facet_grid(parent ~ STIM) +
  stat_compare_means(aes(group = Sex))

In many of the boxplots above, Males are shifted slightly higher on the FS/PFS axis compared to Females, but none of the comparisons are significant. Perhaps a t-test or un-stratified plot would be significant.

Hosp vs Non-Hosp

  1. Plot the COMPASS heatmaps again, this time including row annotations for Cohort
  2. Compare FS and PFS between Hosp and Non-Hosp

COMPASS Heatmaps

crList2 <- lapply(crList,
                  function(cr) {
                    cr$data$meta$Cohort <- factor(ifelse(cr$data$meta$Cohort %in%
                                                c(NA, "Healthy control"), "Healthy", cr$data$meta$Cohort),
                                              levels = rev(c("Healthy", "Non-hospitalized", "Hospitalized")))
                    cr
                  })
names(crList2) <- names(crList)

CohortColors <- c("Healthy" = "#dfc27d", "Non-hospitalized" = "#80cdc1", "Hospitalized" = "#018571")
row_ann_colors <- list(`Cohort` = CohortColors)
cytokine_annotation_colors <- rep("black", 30)
plotCOMPASSResultStack(crList2[CD4RunNames],
                       row_annotation = c("Stim", "Cohort"),
                       variable = "Stim",
                       show_rownames = FALSE,
                       main = "CD4 COMPASS Runs",
                       fontsize = 15, fontsize_row = 15, fontsize_col = 13,
                       row_annotation_colors=row_ann_colors,
                       cytokine_annotation_colors=cytokine_annotation_colors,
                       order_by_max_functionality = T)
## The 'threshold' filter has removed 6 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&IL17a Ax700&!IL2 PE&!IL4/5/13 APC&!TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&!TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&IL4/5/13 APC&!TNFa FITC

plotCOMPASSResultStack(crList2[NotCD4RunNames],
                       row_annotation = c("Stim", "Cohort"),
                       variable = "Stim",
                       show_rownames = FALSE,
                       main = "Not-CD4 COMPASS Runs",
                       fontsize = 15, fontsize_row = 15, fontsize_col = 13,
                       row_annotation_colors=row_ann_colors,
                       cytokine_annotation_colors=cytokine_annotation_colors,
                       order_by_max_functionality = T)
## The 'threshold' filter has removed 2 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&IL17a Ax700&!IL2 PE&IL4/5/13 APC&!TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC

plotCOMPASSResultStack(crList2[c(CD4RunNames, NotCD4RunNames)],
                       row_annotation = c("Stim", "Cohort"),
                       variable = "Stim",
                       show_rownames = FALSE,
                       main = "CD4 and Not-CD4 COMPASS Runs",
                       fontsize = 15, fontsize_row = 15, fontsize_col = 13,
                       row_annotation_colors=row_ann_colors,
                       cytokine_annotation_colors=cytokine_annotation_colors,
                       order_by_max_functionality = T)
## The 'threshold' filter has removed 3 categories:
## !CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&IL2 PE&!IL4/5/13 APC&TNFa FITC, CD107a PE-Cy7&!CD154 PE-Cy5&!IFNg V450&!IL17a Ax700&!IL2 PE&!IL4/5/13 APC&TNFa FITC, !CD107a PE-Cy7&!CD154 PE-Cy5&IFNg V450&!IL17a Ax700&IL2 PE&IL4/5/13 APC&!TNFa FITC

I think a commonality to the subsets which are enriched in Hospitalized individuals is CD154.

Now look at each run individually

cytokine_order_for_annotation = c("IFNg V450", "IL2 PE", "TNFa FITC", "CD154 PE-Cy5", "CD107a PE-Cy7", "IL4/5/13 APC", "IL17a Ax700")

for(cr_name in c(CD4RunNames, NotCD4RunNames)) {
  print(my_plot.COMPASSResult(crList2[[cr_name]],
                        y="Cohort",
                        #show_rownames = T,
                        fontsize=16, fontsize_row=16, fontsize_col=19,
                        row_annotation_colors=row_ann_colors,
                        cytokine_annotation_colors=cytokine_annotation_colors,
                        cytokine_height_multiplier=2.5, order_by_max_functionality = T,
                        main=cr_name,
                        cytokine_order_for_annotation=cytokine_order_for_annotation, 
                        staircase_cytokine_annotation=T # Looks prettier and is consistent with dotplots
                        ))
}
## The 'threshold' filter has removed 7 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.4657]
## The 'threshold' filter has removed 4 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.4785]
## The 'threshold' filter has removed 6 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.4922]
## The 'threshold' filter has removed 10 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.5037]
## The 'threshold' filter has removed 2 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.5164]
## The 'threshold' filter has removed 3 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.5286]
## The 'threshold' filter has removed 3 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7

## gTree[GRID.gTree.5417]
## The 'threshold' filter has removed 6 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7

## gTree[GRID.gTree.5528]

FS/PFS vs Cohort

# Note that p-values on ggboxplot plots below are *not* adjusted, e.g.:
broom::tidy(wilcox.test(FS ~ Cohort, data = fs_pfs_df_with_manifest %>% dplyr::filter(Cohort != "Healthy" & parent == "4" & STIM == "Spike_2")))
## # A tibble: 1 x 4
##   statistic  p.value method                 alternative
##       <dbl>    <dbl> <chr>                  <chr>      
## 1       162 0.000345 Wilcoxon rank sum test two.sided
# Note that future ggpubr update should make outlier.shape=NA occur by default when jitter is added, but include for now
ggboxplot(fs_pfs_df_with_manifest %>% 
            dplyr::filter(Cohort != "Healthy"),
          x = "Cohort", y = "FS",
          color = "Cohort", palette = "jco",
          add = "jitter",
          outlier.shape = NA) + 
  facet_grid(parent ~ STIM) +
  stat_compare_means(aes(group = Cohort)) +
  scale_x_discrete(labels=c("Non-hospitalized" = "Conv\nNon-Hosp", "Hospitalized" = "Conv\nHosp"))

ggboxplot(fs_pfs_df_with_manifest %>% 
            dplyr::filter(Cohort != "Healthy"),
          x = "Cohort", y = "PFS",
          color = "Cohort", palette = "jco",
          add = "jitter",
          outlier.shape = NA) + 
  facet_grid(parent ~ STIM) +
  stat_compare_means(aes(group = Cohort)) +
  scale_x_discrete(labels=c("Non-hospitalized" = "Conv\nNon-Hosp", "Hospitalized" = "Conv\nHosp"))

Background-corrected magnitudes

crList.no_healthy <- lapply(names(crList2),
                  function(cr_name) {
                    cr <- crList2[[cr_name]]
                    print(sprintf("Accessing %s", cr_name))
                    cr$data$meta <- cr$data$meta %>% 
                      dplyr::filter(!(`SAMPLE ID` %in% setdiff(cr$data$meta$`SAMPLE ID`, rownames(cr$fit$mean_gamma))))
                    stopifnot(all.equal(rownames(cr$fit$mean_gamma), cr$data$meta$`SAMPLE ID`))
                    
                    stopifnot(all.equal(rownames(cr$data$n_s), rownames(cr$fit$mean_gamma)))
                    stopifnot(all.equal(rownames(cr$data$n_u), rownames(cr$fit$mean_gamma)))
                    stopifnot(all.equal(names(cr$data$counts_s), rownames(cr$fit$mean_gamma)))
                    stopifnot(all.equal(names(cr$data$counts_u), rownames(cr$fit$mean_gamma)))
                    
                    not_healthy_idx <- which(cr$data$meta$Cohort != "Healthy")
                    cr$data$meta <- cr$data$meta[not_healthy_idx,] %>% 
                      mutate(Cohort = factor(Cohort, levels = c("Non-hospitalized", "Hospitalized")))
                    cr$fit$mean_gamma <- cr$fit$mean_gamma[not_healthy_idx,]

                    cr$data$n_s <- cr$data$n_s[not_healthy_idx,]
                    cr$data$n_u <- cr$data$n_u[not_healthy_idx,]
                    cr$data$counts_s <- cr$data$counts_s[not_healthy_idx]
                    cr$data$counts_u <- cr$data$counts_u[not_healthy_idx]
                    
                    cr
                  })
## [1] "Accessing 4_VEMP"
## [1] "Accessing NOT4_VEMP"
## [1] "Accessing 4_Spike_1"
## [1] "Accessing NOT4_Spike_1"
## [1] "Accessing 4_Spike_2"
## [1] "Accessing NOT4_Spike_2"
## [1] "Accessing 4_NCAP"
## [1] "Accessing NOT4_NCAP"
names(crList.no_healthy) <- names(crList2)

source(here::here("scripts/20200604_Helper_Functions.R"))

# Arial is used by dotplot function
# Arial font setup. Downloaded afms from https://github.com/microsoft/microsoft-r-open/tree/ec3fd89e5fb5794bd8149905c134ad801bb61800
Arial <- Type1Font(family = "Arial",
                   metrics = c(here::here("data/Arial_afm/ArialMT.afm"),
                               here::here("data/Arial_afm/ArialMT-Bold.afm"), 
                               here::here("data/Arial_afm/ArialMT-Italic.afm"),
                               here::here("data/Arial_afm/ArialMT-BoldItalic.afm")))
pdfFonts(Arial = Arial)

# Note p-values are bonferroni adjusted in dotplots
bg_corr_dotplots <- pmap(.l = list(c(CD4RunNames, NotCD4RunNames),
                                   rep(c("CD4+", "Not-CD4+"), each = 4),
                                   list(c(-0.0015, 0.01), c(-0.002, 0.01), c(-0.003, 0.005), NULL,
                                        c(-0.005, 0.007), c(-0.003, 0.008), c(-0.005, 0.013), NULL),
                                   c(5, 5, 3, 5, 5, 5, 5, 5)),
                         .f = function(n, p, y, p_text_size) {
                           # "staircase effect" for categories df heatmap gets done automatically, unlike in my_plot.COMPASSResult
                           # Returned dotplot is cowplot
                           make_dotplot_for_COMPASS_run(crList.no_healthy[[n]], n, parentSubset = p,
                                                        return_output = T, current_ylim = y, include_0_line=T, 
                                                        cytokine_order_for_annotation=cytokine_order_for_annotation,
                                                        p_text_size=p_text_size)
                         })
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
## Joining, by = "BooleanSubset"
names(bg_corr_dotplots) <- c(CD4RunNames, NotCD4RunNames)

for(n in names(bg_corr_dotplots)) {
  print(plot_grid(ggplot() +
    theme_void() +
    geom_text(aes(0,0,label=n), size=10) +
    xlab(NULL),
    bg_corr_dotplots[[n]]$Dotplot,
    ncol=1, rel_heights = c(0.1, 1)))
}

sig_tests_tally <- lapply(names(bg_corr_dotplots), function(n) {length(which(bg_corr_dotplots[[n]]$Test_Results$p.adj < 0.05))})
names(sig_tests_tally) <- names(bg_corr_dotplots)
sig_tests_tally # make sure all significant comparisons were plotted and accounted for
## $`4_Spike_1`
## [1] 8
## 
## $`4_Spike_2`
## [1] 4
## 
## $`4_NCAP`
## [1] 5
## 
## $`4_VEMP`
## [1] 0
## 
## $NOT4_Spike_1
## [1] 0
## 
## $NOT4_Spike_2
## [1] 0
## 
## $NOT4_NCAP
## [1] 0
## 
## $NOT4_VEMP
## [1] 0

Focus on IL2

Look at each run individually again, this time with IL2 subsets on RHS (since there was evidence of elevated DMSO IL2 signal in hosp)

for(cr_name in c(CD4RunNames, NotCD4RunNames)) {
  print(my_plot.COMPASSResult(crList2[[cr_name]],
                        y="Cohort",
                        #show_rownames = T,
                        fontsize=16, fontsize_row=16, fontsize_col=19,
                        row_annotation_colors=row_ann_colors,
                        cytokine_annotation_colors=cytokine_annotation_colors,
                        cytokine_height_multiplier=2.5, order_by_max_functionality = T,
                        cytokine_order_for_annotation = c("IL4/5/13 APC", "IFNg V450", "TNFa FITC", "IL17a Ax700",
                                                          "CD154 PE-Cy5", "CD107a PE-Cy7", "IL2 PE"),
                        dichotomize_by_cytokine = "IL2 PE",
                        staircase_cytokine_annotation = TRUE,
                        main=cr_name))
}
## The 'threshold' filter has removed 7 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7
## gTree[GRID.gTree.8731]
## The 'threshold' filter has removed 4 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.8859]
## The 'threshold' filter has removed 6 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.8996]
## The 'threshold' filter has removed 10 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.9111]
## The 'threshold' filter has removed 2 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.9238]
## The 'threshold' filter has removed 3 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7

## gTree[GRID.gTree.9360]
## The 'threshold' filter has removed 3 categories:
## !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7

## gTree[GRID.gTree.9491]
## The 'threshold' filter has removed 6 categories:
## !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&!IFNg V450&!TNFa FITC&!IL17a Ax700&CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&IL4/5/13 APC&!IFNg V450&!TNFa FITC&IL17a Ax700&!CD154 PE-Cy5&!CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&!TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7, !IL2 PE&!IL4/5/13 APC&IFNg V450&TNFa FITC&!IL17a Ax700&!CD154 PE-Cy5&CD107a PE-Cy7

## gTree[GRID.gTree.9602]

Previously, we observed higher DMSO IL2 signal in Conv Hosp compared to Conv Non-hosp, but so far there’s no obvious sign that this difference carries over universally to the background-corrected IL2-containing subsets. Quantitative magnitude boxplots should provide more detail.